
function [ll,Dll] = LL_cfs(theta,X_PVEM,X_PRI,ES0,ES1,ES2,n)

N0 = n(1,1); N1 = n(1,2); N2 = n(1,3);

P1 = exp(X_PVEM * theta + ES1 - ES0);
P2 = exp(X_PRI * theta + ES2 - ES0);
P0 = 1 + P1 + P2;
P1 = P1 ./ P0;
P2 = P2 ./ P0;
P0 = 1 ./ P0;

ll = sum(log(P0(1:N0))) + sum(log(P1(N0+(1:N1)))) + sum(log(P2(N0+N1+(1:N2))));
ll = -ll;

if nargout > 1
    Dll = sum(([zeros(N0,1);ones(N1,1);zeros(N2,1)]-P1) .* X_PVEM + ...
        ([zeros(N0+N1,1);ones(N2,1)] - P2) .* X_PRI,1)';
    Dll = -Dll;
end



